T    = 1.0d2
y    = 0.0d0
nH   = 9.482359456827762E-4
G0   = 6.317259396184505E-10
tau0 = 220.73506384955874

sigma0 = 6.3d-18

N=1000000

RC = hiirca_hui(T)
CI = hici_hui(T)

R = (CI+RC)*nH
Q = -( (CI+2*RC)*nH + (CI+RC)*nH*y )
P = RC*nH*(1.0d0+y)
S = -G0

; want x (ionized fraction) to go from xHII_0 to small

xHI_min = 1.0d0-6
xHI_max = 1.0d0
xHI_ana = dblarr(N)

xHI_ana = xHI_min + dindgen(N)/(N-1) * (xHI_max - xHI_min) 

lhs = R * xHI_ana*xHI_ana + Q * xHI_ana + P + S * exp(-tau0*xHI_ana) * xHI_ana

altay_set_x
altay_plot, alog10(xHI_ana), lhs

altay_oplot, [-100,100], [0,0], color=100
altay_oplot, [-10,-10], [-1000,1000], color=100










end
